Orbital magnetization in crystalline solids: 
Multi-band insulators, Chern insulators, and metals 
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We derive a multi-band formulation of the orbital magnetization in a normal periodic insulator 
(i.e., one in which the Chern invariant, or in 2d the Chern number, vanishes). Following the approach 
used recently to develop the single-band formalism [T. Thonhauser, D. Ceresoli, D. Vanderbilt, and 
R. Resta, Phys. Rev. Lett. 95, 137205 (2005)], we work in the Wannier representation and find that 
the magnetization is comprised of two contributions, an obvious one associated with the internal 
circulation of bulk-like Wannier functions in the interior and an unexpected one arising from net 
currents carried by Wannier functions near the surface. Unlike the single-band case, where each 
of these contributions is separately gauge-invariant, in the multi-band formulation only the sum of 
both terms is gauge-invariant. Our final expression for the orbital magnetization can be rewritten 
as a bulk property in terms of Bloch functions, making it simple to implement in modern code 
packages. The reciprocal-space expression is evaluated for 2d model systems and the results are 
verified by comparing to the magnetization computed for finite samples cut from the bulk. Finally, 
while our formal proof is limited to normal insulators, we also present a heuristic extension to Chern 
insulators (having nonzero Chern invariant) and to metals. The validity of this extension is again 
tested by comparing to the magnetization of finite samples cut from the bulk for 2d model systems. 
We find excellent agreement, thus providing strong empirical evidence in favor of the validity of the 
heuristic formula. 

PACS numbers: 75.10.-b, 75.10.Lp, 73.20.At, 73.43.-f 



I. INTRODUCTION 

During the last decade, charge and spin transport phe- 
nomena in magnetic materials and nanostructures have 
attracted much interest due to their important role for 
spintronic devices. 1 An adequate description of mag- 
netism in these materials, however, should not only in- 
clude the spin contribution, but also should account 
for effects originating in the orbital magnetization. In 
light of this, it is surprising that the theory of or- 
bital magnetization has long remained underdeveloped. 
Earlier attempts to develop such a theory used linear- 
response methods, which allow calculations of magneti- 
zation changes, 2 ^ but not of the magnetization itself. 

Just recently, a new approach using Wannier functions 
(WFs) has been proposed, 6,7 which nicely parallels the 
analogous case of the electric polarization. The primary 
difficulty in both cases is that the position operator r is 
not well-defined in the Bloch representation. Since WFs 
are exponentially localized in an insulator, this difficulty 
disappears if the problem is reformulated in the Wannier 
representation. For the polarization, this approach lead 
to the development of the modern theory of polarization 
in the early 1990s. 8 ' 9 Similarly, in the case of the orbital 
magnetization, where the circulation operator r x v is 
ill-defined in the Bloch representation, the Wannier rep- 
resentation was used to derive a theory for the orbital 
magnetization of periodic insulators. 7 

While the formalism developed in Ref. 7 lays a firm 
foundation for the orbital magnetization, its application 
is limited to certain systems, such as single-band models 



and insulators. In this paper we expand the applica- 
bility to a much wider class of systems by developing a 
corresponding multi-band formalism, essential for most 
"real" materials. This extension is nontrivial and the 
corresponding proof of gauge invariance is much more 
complex than for the single-band case. We proceed in 
two steps. First, we carry out a derivation for the case of 
an insulator with zero Chern invariant. Second, we give 
heuristic arguments for an extension of our formalism to 
metals and Chern insulators, i.e. systems with a non-zero 
Chern invariant, arriving at a formula identical to that 
proposed by Xiao, Shi and Niu 10 on the basis of semiclas- 
sical arguments. Chern insulators have been introduced 
into the theoretical literature by means of model Hamil- 
tonians in 2d which break time-reversal (TR) symmetry 
without breaking translational symmetry, 11 i.e., main- 
taining a vanishing macroscopic magnetic field. Despite 
the absence of a macroscopic field, Chern insulators share 
several properties with quantum-Hall systems, most no- 
tably the quantization of the transverse conductivity in 
2d. 11 To the best of our knowledge, there is no known ex- 
perimental realization of a Chern insulator (in zero field) 
in either 2d or 3d, and the search for such a system re- 
mains a fascinating challenge. 

Our extensions to metals and Chern insulators are 
heuristic and not based on an analytical proof. The fact 
that our final formula is identical to the one derived from 
the semiclassical wavepacket treatment 10 is reassuring, 
but neither of these approaches can yet be said to consti- 
tute a "derivation" of the formula in the fully quantum 
context. Nevertheless, we provide strong numerical evi- 
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dence of their validity, thus posing a theoretical challenge: 
how to provide an analytic proof of the heuristic formula, 
beyond the range of the semiclassical approximation, for 
both the metallic and Chern-insulating cases. 

Before proceeding, we emphasize that the present work 
only addresses the question of how to compute the orbital 
magnetization for a given independent-particle Hamil- 
tonian. Many interesting questions remain concern- 
ing which flavor of density-functional theory (DFT) or 
which exchange-correlation (XC) functional might give 
the most accurate orbital magnetization. While exact 
Kohn-Sham (KS) density (or spin-density) functional 
theory is guaranteed to yield the correct charge (or spin) 
density, 12 there is no reason to expect it to yield the 
correct orbital currents. The orbital magnetization, be- 
ing defined in terms of surface currents, is not guar- 
anteed to be correct either. A prescription that seems 
more suited to the present situation is that of Vignalc 
and Rasolt, 13 in which the spin-labeled density and cur- 
rent {n CT (r), jcr(r)} are connected to corresponding scalar 
and vector potentials {V a (r), A cr (r)}. However, it is an 
open question whether an approximate Vignale-Rasolt 
XC functional exists that can give improved values of 
magnetization in practice. While time-dependent den- 
sity functional theory (TDDFT) is more developed, 14 
this theory only establishes a connection between n(r, t) 
and V(r,t), and a knowledge of n(r, t) is only suffi- 
cient to determine the longitudinal part of j(r, t), not 
the transverse part upon which the orbital magnetiza- 
tion depends. An alternative approach worthy of explo- 
ration is time-dependent current-density functional the- 
ory (TD-CDFT), 15 in which {n(r, t), j(r, t)} is connected 
to {V(r, t), A(r, t)}. However, the present problem is 
essentially a static problem, and it is therefore unclear 
whether TD-CDFT would provide any practical advan- 
tage over the Vignale-Rasolt theory. Finally, it is worth 
remembering that even in standard DFT, the mapping 
from interacting density to non-interacting potential is 
sometimes pathological (e.g., a KS metal can represent 
an interacting insulator). In the present work, we bypass 
all these interesting issues, and only consider how to com- 
pute the magnetization for a given Kohn-Sham Hamilto- 
nian arising from some unspecified version of DFT in the 
context of broken TR symmetry. 

We have organized this paper as follows. In Sec. II we 
derive the multi-band theory of orbital magnetization in 
crystalline solids. After some definitions and generalities, 
we start by considering the orbital magnetization of a fi- 
nite sample. The resulting expression is then transformed 
to reciprocal space and its gauge invariance is demon- 
strated. We then give a heuristic extension of our formal- 
ism to metals and Chern insulators. In Sec. Ill, numeri- 
cal results for the orbital magnetization are presented for 
several different systems. We conclude in Sec. IV. Some 
details concerning the finite-difference evaluation of the 
magnetization and certain properties of the nonAbelian 
Berry curvature are deferred to two appendices. 



II. THEORY 
A. Generalities 

Our basic starting point is a single-particle KS Hamil- 
tonian 12 having the translational symmetry of the crys- 
tal, but having no TR symmetry: as said above, trans- 
lational symmetry of the Hamiltonian implies vanishing 
of the macroscopic magnetic field. There may, however, 
be a microscopic magnetic field B that averages to zero 
over the unit cell, and we assume that a particular mag- 
netic gauge has been chosen once and for all to represent 
this magnetic field. Wavevector k is a good quantum 
number under these conditions. This could be realized, 
for example, in systems in which the TR breaking comes 
about through the spontaneous development of ferromag- 
netic order or via spin-orbit coupling to a background of 
ordered local moments. 11,16-19 Notice that we carefully 
avoid referring to an externally applied field; such con- 
cept is legitimate only for a finite sample, free-standing 
in vacuo. Indeed, for a finite sample, the relationship 
between the externally applied field and the "internal" 
(or screened) one depends on the sample shape. For an 
extended sample in the thermodynamic limit, the only 
legitimate and measurable field is the screened B field 
which is present inside the material. In the present work, 
the cell-average of this field is assumed to vanish. 

As usual, we let e„k and \ip n k) be the Bloch eigenval- 
ues and eigenvectors of H, respectively, and w„k(r) = 
C k ' r ^nk(i') be the corresponding eigenfunctions of the 
effective Hamiltonian = e _lk r iJe lk r . We choose to 
normalize them to one over the crystal cell of volume Q. 

The notation is intended to be flexible as regards the 
spin character of the electrons. If we deal with spinless 
electrons, then n is a simple index labeling the occupied 
Bloch states; factors of two may trivially be inserted if 
one has in mind degenerate, independent spin channels. 
In the context of the local spin-density approximation 
(LSD A), in which spin-up and spin-down electronic states 
are separate eigenstates of spin-up and spin-down Hamil- 
tonians, one may let n range over both sets of bands, but 
with the understanding that inner products or matrix ele- 
ments between spin-up and spin-down bands always van- 
ish. Of more realistic interest here is the case of a fully 
non-collinear treatment of the magnetism, as for the case 
of a Hamiltonian containing the spin-orbit operator. In 
this case, n labels bands that are neither purely spin-up 
nor spin-down, |u„k) must be understood to be a spinor 
wavefunction, and the contraction over spin degrees of 
freedom is understood to be included in the definition of 
inner products like (Mnk|w n 'k) and matrix elements like 
(u nk \H k \u n ' k ). 

A key issue in the present work is the additional "gauge 
freedom" in which the occupied Bloch orbitals at fixed 
k are allowed to be transformed among themselves by 
an arbitrary unitary transformation. In fact, any KS 
ground-state electronic property should be uniquely de- 
termined by the subspace of occupied orbitals as repre- 
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sented by the one-particle density matrix; the occupied 
orbitals just provide a convenient orthonormal represen- 
tation for this subspace. Moreover, when it comes to the 
formulation of Wannicr functions (WFs) for composite 
energy bands, the n-th WF is generally not simply the 
Fourier transform of the n-th band of Hamiltonian eigen- 
vectors, but instead, of a manifold of states |u„k) which 
are related to the eigenstates by a k-dependent unitary 
transformation. 20 Thus, in what follows, we allow |u n k) 
to refer to this generalized interpretation of the nk la- 
bels unless otherwise specified. In addition, we introduce 
a generalized "energy matrix" 



E nn 'k — (u n k\H k \u n >k). 



(1) 



which reduces to E nn >k = e n k8 nn i in the special case of 
the "Hamiltonian gauge" in which the |u„k) are eigen- 
states of Hk. 

A key quantity characterizing a three-dimensional KS 
insulator in absence of TR symmetry is the (vector) 
Chern invariant 21 



C = ^— dk V7<9 k Unk| x |<9k«„k), 
2n Jbz „ 



(2) 



with the usual meaning of the cross product between 
three-component bra and ket states. Here and in the 
following the sum is over the occupied n's only, the in- 
tegral is over the Brillouin zone (BZ), and 9k = d/dk. 
The Chern invariant is gauge-invariant in the above gen- 
eralized sense (as will be shown in Sec. II D) and — for 
a three-dimensional crystalline system — is quantized in 
units of reciprocal-lattice vectors G. In Sees. IIB-IID 
we assume that we are working with insulators with zero 
Chern invariant; the more general case will be discussed 
only later in Sees. IIE-IIF. 

Owing to the zero-Chern-invariant condition, the 
Bloch orbitals can be chosen so as to obey |i/Vik+G) = 
iVVik) (the so-called periodic gauge), which in turn war- 
rants the existence of Wannier functions (WFs) enjoying 
the usual properties. (For a Chern insulator, it is not 
clear whether a Wannier representation exists.) We shall 
denote as |nR) the n'th WF in cell R. These WFs are 
related via 

Kk} - 5> ik -( R - r )|nR), 

R 

|nR) = f dkc^-^Kk), (3) 

to the Bloch-likc orbitals |u„k) defined in the generalized 
sense discussed just above Eq. (1). 



B. The magnetization of a finite sample 

We start by considering a macroscopic sample of N c 
cells (with iV c very large but finite) cut from a bulk insu- 
lator, having Nb occupied bands, with "open" boundary 



conditions. The finite system then has N ~ N c Nb occu- 
pied KS orbitals. Suppose we perform a unitary transfor- 
mation upon them, by adopting some localization crite- 
rion. By invariance of the trace the orbital magnetization 
of the finite system is written in terms of the localized 
orbitals \wi) as 



N 



M 



2cQ,N, 



(4) 



i=i 



where the velocity is defined as 

v = i[H,r}. (5) 

In the case of density-functional implementations, it 
should be noted that v may differ from p/m because of 
the presence of microscopic magnetic fields (which intro- 
duce p • A terms in the Hamiltonian), spin-orbit interac- 
tions, or semilocal or nonlocal pseudopotentials. In the 
case of tight-binding implementations, the matrix rep- 
resentations of H and r are assumed to be known (r is 
normally taken to be diagonal) in the tight-binding basis, 
and v is then defined through Eq. (5). 

We divide the sample into an "interior" and a "sur- 
face" region, in such a way that the latter occupies a 
non-extensive fraction of the total sample volume in the 
thermodynamic limit. The orbitals |tOj) which are local- 
ized in the interior region converge exponentially to the 
WFs |nR) of the periodic infinite system; for instance, if 
the Boys 22 localization criterion is adopted, they become 
by construction the Marzari-Vandcrbilt 20 maximally lo- 
calized WFs. Therefore the interior is composed of an 
integer number Ni of replicas of a unit cell containing 
Nb WFs each. Note that this choice is not unique; there 
is freedom both to shift all of the R's by some constant 
vector (effectively changing the origin of the unit cell), 
or to shift any one of the WFs by a lattice vector, or to 
carry out a unitary remixing of the bands. We insist only 
that some consistent choice is made once and for all. 

The remaining N s localized orbitals residing in the sur- 
face region need not resemble bulk WFs; we denote them 
as \w s ) and continue to refer to them as "WFs" in a gen- 
eralized sense. We thus partition the entire set of N WFs 
of the finite sample into NiNb ones belonging to the inte- 
rior and N s ones in the surface region. Correspondingly, 
the contribution to the orbital magnetization M coming 
from the interior orbitals will be denoted as Mlc (for 
"local circulation"), while that arising from the surface 
orbitals will be referred to as Mic (for "itinerant circu- 
lation"). We will take the thermodynamic limit in such 
a way that N s grows more slowly with sample size than 
does Ni, so that N s /Ni — ► 0. Because of the ambiguities 
discussed in the previous paragraph, we do not expect 
Mlc and Mic to be separately gauge-invariant. How- 
ever, their sum, Eq. (4), must be gauge- invariant. 

Since the interior orbitals are bulk-like, we have, fol- 
lowing Eq. (4), 

M L c = - ;tt7777 X^ nR K r ~ R ) X v l nR )i ( 6 ) 

nR 



2cVLN, 



4 



s x s 

FIG. 1: Horizontal slice from a sample that extends indefi- 
nitely in the vertical direction. Vertical dashed lines delimit 
bulk and surface regions in which WFs are labeled by s and 
s' , respectively. 



where the number of R vectors in the sum is smaller than 
N c only by a nonextensive fraction, and we have used 
that ^ n (nR|v|nR) = 0. Because of the zero-Chern- 
invariant condition the WFs enjoy the usual translational 
symmetry, and we finally find that 



M LC = - ^(nO|r x v|nO) 



(7) 



in the thermodynamic limit. 

We now consider the contribution from the N s surface 
orbitals, whose centers we denote as r s — (w s \r\w s ) : 



Mic 



1 



2cflN, 



E 



((w s \(r-r s )xv\w s )+r s x(w s \\r\w s )). 

(8) 

The first term in parenthesis clearly vanishes in the ther- 
modynamic limit, while the second term — owing to the 
presence of the "absolute" coordinate r s — does not. At 
first sight, this second term in Mjc appears to depend on 
surface details; instead, we are going to prove that even 
this term can be recast in terms of bulk Wannier func- 
tions. Remarkably, both Mlc and Mic are genuine bulk 
properties in the thermodynamic limit, and can eventu- 
ally be evaluated as BZ integrals. 

We consider a surface facing in the +x direction, and 
identify a surface region given by x > xq as in Fig. 1. 
There is then a contribution to the macroscopic surface 
current K flowing at the surface that is given by 



|v|uv), 



(9) 



where the primed sum is taken over the surface WFs 
whose yz coordinates are within one surface unit cell of 
area A. Because (uv|v|uv) decays exponentially to zero 
with distance from the surface, it is straightforward to 
capture the entire surface current by letting the width of 
the surface region diverge slowly (say, as the 1/4 power 
of linear dimension) in the thermodynamic limit, so that 
xo is moved arbitrarily deep into the bulk. 
It is now expedient to use the identity 



(Wi\v\wi) = J^v 



where 



'CM) 



21m (wi\r\wj)(wj\H\wi) 



(10) 



(11) 



has the interpretation of a current "donated from WF 
\wj) to WF \wi)" , and exploit the fact that the total 
current carried by any subset of WFs can be computed 
as the sum of all v<j t i\ for which i is inside and j is outside 
the subset. Applying this to the piece of surface region 
considered above, we get 



(12) 



Setting the boundary deep enough below the surface to 
be in a bulk-like region and invoking the exponential lo- 
calization of the WFs and of related matrix elements, we 
can identify \w s ) and \w S '} with the bulk WFs |mR) and 
|nR'), respectively. Exploiting translational symmetry, 
V( m R,„R'> = v (m0 ,„(R/_R)), Eq. (12) becomes 

1 X - X - X - 



K = — - 



A ^ ^ ^ 

R x <x R' >x mn 



V( m o,n(R'-R)) 



(13) 



where the lattice sum is still restricted to the R' vectors 
whose yz coordinates are within the surface unit cell. The 
number of terms in the lattice sum of Eq. (13) having a 
given value ofR'-Risjust (R' X -R X )A/Q, if (R' x -R x ) > 
and zero otherwise. With a change of summation index, 
Eq. (13) becomes 



R mn 



v (mO,nR) j 



(14) 



where the factor of 2 enters because the sum has been 
extended to all R. Notice that the surface-cell size has 
eventually disappeared. 

Evidently the corresponding surface current on a sur- 
face with unit normal n is then 



K a (n) 







where 



G, 



a/3 



R mn 



"(mO,nR),a-ft/3 



(15) 



(16) 



Now for a sample of size L x x L v x L z , the left and 
right faces carry currents of ±L y L z G yx separated by a 
distance L x , and thus contribute to the magnetic moment 
per unit volume as G yx /2c; similarly, the front and back 
faces contribute as —G xy /2c. Together they contribute 
to M z as —G x /c where 



Gap — 2 ( ^"z 3 _ Gj3a 



(17) 



is the antisymmetric part of the G tensor. Deriving cor- 
responding expressions for M x and M y by permutation of 
indices, the contribution of the surface current in Eq. (14) 
to the magnetization can thus be cast in a coordinate- 
independent form and evaluated for the whole sample 
surface in the thermodynamic limit as 



Mic 



1 



X V 



(mO,nR) 



(18) 



mnR 



■5 



Note that Eq. (18) describes the current circulating in 
the surface WFs, while the expression on its right-hand 
side involves only bulk WFs. 

This is quite remarkable, and indeed it is one of the 
central results of this paper, as well as of Ref. 7. It im- 
plies that even Mic is a bulk property, as anticipated 
above. This may appear counterintuitive, but indeed 
closely parallels a well-known (and equally counterintu- 
itive) feature of the quantum-Hall effect, where the Hall 
current is accomodated by chiral edge states. 23,24 Never- 
theless, these edge currents are completely determined by 
bulk properties of the system, and can be evaluated by 
adopting toroidal boundary conditions in which the sam- 
ple has no edges. Such a finding, in fact, is one of the 
most remarkable results of the quantum-Hall theoretical 
literature. 21,25-27 We also notice that the bulk nature of 
Mic guarantees that our general expressions, valid in the 
thermodynamic limit, apply regardless of whether sur- 
face states are present in bounded samples, and if they 
are present, regardless of their character. 

It might be thought that the surface currents K must 
flow parallel to the surface, and thus that the diagonal 
elements G xx and G yy must vanish, or more generally, 
that the symmetric component 



G 



af3 — 9 ( G »/3 + Gpa 



(19) 



of the G tensor should vanish. This turns out not to be 
true. In some of our tight-binding model calculations, we 
have explicitly computed the right-hand side of Eq. (9) 
and confirmed the existence of a surface-normal compo- 
nent of K. 

The explanation is that K, as defined by Eq. (9), is 
only one contribution to the physical macroscopic sur- 
face current. There is an additional contribution arising 
from the fact that, when TR symmetry is broken, the 
second-moment spreads 20 of the WFs are not generally 
stationary with respect to time. For example, if the WFs 
are in the process of expanding, then electron charge is 
in the process of spilling out of the surface. To formalize 
this notion, we introduce a symmetric Cartesian tensor 



W, 



1 



a/3 



2flN r ^ 



r a vp + v a r fi \w l ) (20) 



that is a kind of symmetric analog of the antisymmet- 
ric expression for M given in Eq. (4). If W a fj is non- 
zero, then we would expect surface currents of the form 
K a (n) = Waprifj. If present, these would violate con- 
tinuity. However, they are not present, because we can 
write 



W af} = - 



2flN r dt 



^2(wi\r a rp\wi) . (21) 



Noting that the trace of any operator (here r a rp) must 
be independent of time in any stationary state (here the 
ground state of the finite sample) , it follows that W a p — 



0. Nevertheless, if we were to follow a route parallel to 
that used for the treatment of M earlier in this section, 
we could decompose W into a "local spread" part Wls 
and an "itinerant spread" part Wis- The former is 

W LS , a p = --^r^2(n0\r a v p +v a rf3\n0) 



2Q 
1 d 



(22) 



which is just related to the rate of spread of the bulk WFs 
in one bulk unit cell, while the latter is just Wis, a/3 = 
G^p of Eq. (19). Because the total W a p must vanish, 
we conclude that the non-physical current that we were 
concerned about, arising from G| in Eq. (19), is exacly 
cancelled by another non-physical one arising from the 
spreading of the bulk WFs. Thus, in the end, the physical 
edge current has pure circulating character and is related 
only to antisymmetric Cartesian tensors. 



C. Reciprocal-space expressions 

The above two final expressions, Eqs. (7) and (18), 
are given in terms of bulk WFs. Therefore the total 
orbital magnetization M = Mlc + Mic of the finite 
sample converges in the thermodynamic limit to a bulk, 
boundary- insensitive, material property. Next, using the 
WF definition, Eq. (3), we are going to transform Mlc 
and Mic into equivalent expressions involving BZ inte- 
grals of Bloch orbitals. Specifically, we are going to prove 
the two identities 



M LC 



1 



2c(2tt) ; 



IrnE 



dk (<9kU Ilk | x iJ k |<9 k u„ k ) 

BZ 



(23) 



U n 'k) ■ 

(24) 

These two expressions generalize to the multi-band case 
our previous finding for the case of a single occupied 
band. 7 There is an important difference, however; while 
in the single-band case Eqs. (23) and (24) are separately 
gauge- invariant, only their sum is gauge-invariant in the 
multi-band case, as we shall see in Sec. II D. 

We carry the derivation in reverse, starting from 
Eqs. (23) and (24) and showing that they reduce to 
Eqs. (7) and (18). First, using Eq. (3), we get 

|cWk) = -^c jk '( R - r )(r-R)|nR) , 

R 

Hk|cW k ) = -^c lk '( R - r )ff(r-R)|nR) . (25) 

R 

Since the velocity operator is v = i[H, r] = i[H, (r — R)], 
and exploiting (r — R) x (r — R) = 0, we may express 
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Eq. (23) as 

MLC = ~2d^^ (nR|(r ~ R)Xv|nR> ' (26) 



TlR 



where the number of cell iV c here is formally infinite, 
and appears because the \u n k) are normalized differently 
from the WFs. Since we limit ourselves to the case of 
an insulator with zero Chern invariant, the WFs enjoy 
the usual translational symmetry, and Eq. (26) is indeed 
identical to Eq. (7). 

Next, we address Eq. (24), whose second factor in the 
integral is 

(<9kWnk| X |<9kU„' k ) = 

= -i- e lk ' (R '~ R) (nR| (r - R) x (r - R') |n'R') 



RR' 



1 J2 gik.(R'-R) (R , _ R) x ^ R | r | n / R /^ (27) 



RR' 



where the last line follows because only the cross terms 
survive from the product (r — R) x (r — R'). We then 
exploit 



(n'R'\H\nR) = / dk e* k < R '-^E n , n k 



(28) 



in order to rewrite Eq. (24) as 
Im 



M 



ic 



2cClN c 



^(R'-R) x (raR|r|nR')(nR'|i2JmR). 



RR' 

ran 



(29) 

Since the matrix elements only depend on the relative 
WF coordinate R' — R, Eq. (29) is transformed into 

Mic = -^r Im ^Rx (™0|r|nR)(nR| J ff|mO). (30) 



mtiR 



Using Eq. (11), it is then easy to check that Eq. (30) is 
indeed identical to Eq. (18). 

This completes our proof. Our final expression for the 
macroscopic orbital magnetization of a crystalline insu- 
lator is 



M 



x ( ff k <W + E n , nk ) |c\u nk ) . (31) 



Owing to the occurrence of Hk and £wk with the same 
sign (in contrast to the magnetization of an individual 
wavepacket discussed in Ref. 28), Eq. (31) does not ap- 
pear at first sight to be invariant with respect to trans- 
lation of the energy zero. However, the zero-Chcrn- 
invariant condition — compare Eq. (31) to Eq. (2) — 
enforces such invariance. As for the gauge invariance of 
Eq. (31), this will be demonstrated in the next subsec- 
tion. 



D. Proof of gauge invariance 

Here we prove the gauge invariance in the multi-band 
sense of the Chern invariant, Eq. (2), and of our main 
expression for the macroscopic magnetization, Eq. (31). 
While these expressions are BZ integrals, we will actually 
prove that even their integrands are gauge-invariant. To 
this end, we will show that both integrands can be ex- 
pressed as traces of gauge-invariant one-body operators 
acting on the Hilbert space of lattice-periodical functions. 

Our key ingredients are the effective Hamiltonian Hk, 
the ground-state projector 



Pk = ^2 \u n u)(unk\ 



(32) 



and its orthogonal complement Qk = 1 — Pk- These 
three operators are obviously unaffected by any unitary 
mixing of the |u„k) among themselves at a given k, and 
therefore any expression built only from these ingredients 
will be a manifestly multi-band gauge-invariant quantity. 
In particular, we define the three quantities 



/k,a/3 = tr{ (d a Pk)Qk (cVk)} 



(33) 



3k,a/3 = tr{ (d a Pk) Qk Hk Qk {dpPk) } , (34) 



Ik, a/3 



tr{tf k (d a Pk) Qk (dpPk)} 



(35) 



where d a — d/dk a and the trace is over electronic states. 
We are going to show that the Chern invariant and the 
magnetization can be expressed as integrals of fk and of 
.9k + ^k, respectively. 

From Eq. (32) it follows that 



d a Pk = ^ ( \ d aU n k)(Unk\ + \u n k) {9 a U n k\ 



(36) 



so that 

{d a Pk)Qk (dpPk) = ^ l U nk)(9 Q U„k|Qk|9/3M„'k)(Wn'k|- 

nn' 

(37) 

We now specialize to the "Hamiltonian gauge" in which 
the Bloch functions |u„k) are eigenstates of Hk with 
eigenvalues e„k- Inserting Eq. (37) into Eqs. (33) and 
(35) and using a similar approach for Eq. (34), the three 
quantities can be written as 

fk.ap = y^jdgUnkldpUnk) 

n 

- ^2(d a u n k\u n 'k)(u n 'k\d/3U n k), (38) 

nn' 

9k, a p = ^2(d a u n k\Hk\d/3U n k) 

n 

- ^2en>k(d a u n k\u nl k)(u n >k\df3U n k), (39) 



7 



and 

hu,a/3 = ^ enk(daU n u\dpU n k) 

n 

~ enk(9 Q W„k|Wn'k)(Un'k|9 / 3U n k)- (40) 

nn' 

Regarded as 3x3 Cartesian matrices, Eqs. (33-35) are 
clearly Hermitian, so that the antisymmetric parts of 
Eqs. (38-40) are all pure imaginary. Thus, the informa- 
tion content of the antisymmetric part of /k, a /3 is con- 
tained in the gauge-invariant real vector quantity 

/k,a = -Inr£ Q( 3 7 /k,/3 7 , (41) 

where e Q/ 3 7 is the antisymmetric tensor. We define g^ a 
and /ik,a hi the corresponding way in terms of gk,/3j 
and /ik,/3 7 respectively. Looking at the second term in 
Eq. (38) and using 9 a (unk|un'k) = d a S nn ' — 0, we find 
that its antisymmetric part vanishes, and in fact fk is 
nothing other than the Berry curvature. We thus recover 
the Chern invariant of Eq. (2) in the form 

C = -L / dk f k . (42) 
27T Jbz 

Next, inspecting the second terms of Eqs. (39) and 
(40), we find that neither of these terms vanishes by it- 
self under antisymmctrization. However, the sum of these 
two terms does vanish under antisymmetrization. Us- 
ing the sum only, and comparing with Eq. (31), we find 
that the magnetization may be written in the manifestly 
gauge-invariant form 

M = M>krh ^ + ^- (43) 

(The sign reflects the fact that the electron has negative 
charge.) This completes the proof that the integrand in 
Eq. (31) is multi-band gauge- invariant. 

Notice that if we take the first term only in Eq. (39) 
and antisymmctrize, we get the integrand in Mlc (times 
a multiplicative constant); the same holds for Eq. (40) 
and Mio However, the second terms in Eqs. (39) and 
(40) have nonzero antisymmetric parts which are essen- 
tial to their gauge-invariance. Therefore, Mlc and Mic 
as defined above are not separately gauge-invariant, ex- 
cept in the single-band case. 7 

On the other hand, it is possible to regroup terms and 
write M = M LC + Mi C , where 

™- = ^/ B f* k (44) 

and 

™ lc = M^rLt k ^ (45) 

are individually gauge- invariant, even in the multi-band 
case. This raises the fascinating question as to whether 



these two contributions to the orbital magnetization are, 
in principle, independently measurable. On the one 
hand, Berry has emphasized in his milestone paper 29 
that any gauge-invariant property should be potentially 
observable. On the other hand, any measurement of or- 
bital magnetization — or, equivalently, of dissipationless 
macroscopic surface currents — will only be sensitive to 
their sum. At the present time we have no insight into 
how to propose an experiment that could distinguish 
them, and we therefore leave this as an open question. 

In Appendix A, we show how to compute fk, gk, and 
hk on a 3D k-mesh using finite-difference methods to 
approximate the derivatives in Eqs. (33-35). 



E. Heuristic extension to metals and Chern 
insulators 

All of the above results are derived under the hy- 
pothesis that the crystalline system is a KS insulator in 
which the Chern invariant, Eq. (2), is zero. These con- 
ditions, in fact, are essential for expressing any ground- 
state property in terms of WFs. Nonetheless the inte- 
grand in our final reciprocal-space expression, Eq. (31), is 
gauge- invariant. This suggests a possible generalization 
to Chern insulators (defined as insulators with nonzero 
Chern invariant) and even to KS metals. 

We notice that Eq. (31) is somehow reminiscent of the 
Berry-phase formula appearing in the modern theory of 
electrical polarization. 8 ' 9 There is an important differ- 
ence, however. In the electrical case, the integrand is 
not gauge-invariant, and the formula corresponding to 
our Eq. (31) only makes sense when integrated over the 
whole BZ, i.e., for a KS insulator. Indeed, macroscopic 
polarization is a well-defined bulk property only for in- 
sulating materials. 30 Instead, orbital magnetization is a 
phenomcnologically well-defined bulk property for both 
insulating and metallic materials. Therefore, it is worth- 
while to investigate heuristically the validity of an ex- 
tension of Eq. (31) to the metallic case, even though we 
cannot yet provide any formal proof. Additionally, we 
also heuristically investigate Chern insulators. Metals 
and Chern insulators share the property that their mag- 
netization has a nontrivial dependence on the chemical 
potential ji. 

We already observed that Eq. (31) is invariant by 
translation of the energy zero, but this owes to the facts 
that the integration therein is performed over the whole 
BZ, and that the Chern invariant is zero. If we aban- 
don either of these conditions, the formula has to be 
modified in order to restore the invariance. To this 
end, we first need to restrict our formulation to the 
"Hamiltonian gauge" , where the energy matrix is diag- 
onal: £„„'k The |u„k) are therefore eigen- 
states of ifk, and the only gauge freedom allowed is now 
the arbitrary choice of their phase. 

In the general case, including metals and Chern insu- 
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lators, we propose to generalize Eq. (31) to 
M = 



2c(2 



l~— ImV / dk (dkU nk \ 



x{H u + e„ k - 2/i) |5kM„ k ) , (46) 



where /z is the chemical potential (Fermi energy). 
Eq. (46) has the desirable invariance property addressed 
above. Furthermore, in the metallic case, Eq. (46) pro- 
vides a magnetization dependent on /x, as it should. In 
the insulating case, when /x is varied in the gap, M 
changes linearly only if the Chern invariant is nonzero, 
and remains constant otherwise. In fact, Eqs. (2) and 
(46) imply that 



dM 

d\i 



1 



c(2tt) 2 



(47) 



for any insulator and \x in the gap. 

The modification from Eq. (31) to Eq. (46) is the min- 
imal one enjoying the desired properties. Furthermore, 
in the single-band case it is essentially identical to a for- 
mula recently proposed by Niu and coworkers, 10 whose 
derivation rests upon semiclassical wavepacket dynamics. 
We provide strong numerical evidence that this formula 
retains its validity well beyond the semiclassical regime, 
and is in fact the exact quantum-mechanical expression 
for the orbital magnetization (in a vanishing macroscopic 
B field). 

An expression related — though not identical — to 
Eq. (46) occurs in the theory of the Hall effect. Upon 
replacement of the quantity in parenthesis with the iden- 
tity, one obtains something proportional to the integral 
of the Berry curvature over occupied portions of the BZ. 
This quantity corresponds to the entire Hall conduc- 
tivity in quantum-Hall systems 25 ' 26 (which are in fact 
two-dimensional Chern insulators 31 ) and the so-called 
"anomalous" Hall term in metals with broken TR sym- 
metry. The theory of the anomalous Hall effect has at- 
tracted much attention in the recent literature. 17 ' 19 ' 32 



cut from a Chern insulator. Owing to the main equation 
V x M = j/c, a macroscopic current of intensity / = cM 
circulates at the edge of any two-dimensional uniformly 
magnetized sample, hence Eq. (49) yields 



dl 
d\x 



C_ 
'2tt' 



(50) 



This is just what is to be expected: raising the chemi- 
cal potential by d[i fills dk/2ir states per unit length, i.e., 
dl = — vdk/2ir; but the group velocity is just v = d/i/dk. 
Thus, Eq. (50) follows with the interpretation that C is 
the excess number of chiral edge channels of positive cir- 
culation over those with negative circulation. Remark- 
ably, the above equations state that the contribution of 
edge states is indeed a bulk quantity, and can be eva- 
luted in the thermodynamic limit by adopting periodic 
boundary conditions where the system has no edges. As 
already observed, this feature may look counterintutitive, 
but a similar behavior has been known for more than 20 
years in the theory of the quantum-Hall effect. 23-26 

In contrast to our case, a magnetic field is usually 
present in the standard theory of the quantum-Hall ef- 
fect, although it is not strictly needed. 11 The role of chi- 
ral edge states is elucidated, for example, 23 ' 24 by con- 
sidering a vertical strip of width £, where the currents 
at the right and left boundaries are ±1. The net cur- 
rent vanishes insofar as fi is constant throughout the 
sample. When an electric field £ is applied across the 
sample, the right and left chemical potentials differ by 
Ait = El and the two edge currents no longer cancel. 
Our Eq. (50) is consistent with the known quantum-Hall 
results. In fact, according to Eq. (50), the net current is 
AI ~ —CA[i/2ir, while the transverse conductivity is de- 
fined by AI = <jt££. We thus arrive at ctt = —C/2n (or, 
in ordinary units, ctt = —Ce 2 /h), which is indeed a cele- 
brated result. 21 ' 25 ~ 27 We stress that the Chern number C 
is a bulk property of the system, and can be evaluated by 
adopting toroidal boundary conditions, where the edges 
appear to play no role. 



F. The two-dimensional case 

In two dimensions, the magnetization is a pseudoscalar 
M, and the Chern invariant is the Chern number C (a 
dimensionless integer) 21 . Our heuristic formula, Eq. (46), 
then becomes 

x(# k +enk-2n)\d k u nk ) . (48) 
The two-dimensional analogue of Eq. (47) is 

^ = -^ (49) 
dfi 2^c' 1 ' 

The physical interpretation of this equation is best under- 
stood in terms of the chiral edge states of a finite sample 



III. NUMERICAL TESTS 

In a previous paper 7 we tested Eq. (48) for the in- 
sulating C = single-band case on the Haldane model 
Hamiltonian, 11 described below (Sec. IIIC). In this spe- 
cial case, Eq. (48) is not heuristic, since we provided an 
analytical proof. We addressed finite-size realizations of 
the Haldane model, cut from the bulk; our analysis con- 
firmed that Mlc arises entirely from the magnetization 
of bulk WFs in the thermodynamic limit, whereas M\c 
arises from current-carrying surface WFs. Both terms 
have also been evaluated in terms of bulk Bloch orbitals, 
by means of Eq. (48), confirming that the orbital mag- 
netization is indeed a genuine bulk quantity. 

Here we extend this program of checking the correct- 
ness of our analytic formulas by carrying out numerical 
tests on our new multi-band formula, Eq. (31), derived 



for the C — insulating case. This is done using a four- 
band model Hamiltonian on a square lattice as described 
below (Sec. Ill A). Furthermore, we perform computer 
experiments to assess whether our hypothetical Eq. (48), 
proposed to cover also the metallic and the C ^ insu- 
lating cases, is consistent with calculations on finite sam- 
ples. We do this for metals in Sec. IIIB using the same 
square lattice as in Sec. Ill A, but at fractional band fill- 
ing. We then do this in Sec. Ill C for Chern insulators 
using the Haldanc model 11 in a range of parameters for 
which C ^ 0. 

Numerical implementation of Eqs. (31), (46), and (48) 
is quite straightforward once one has in hand an efficient 
method for evaluating the k-dcrivatives of the Bloch or- 
bitals. There are several possible approaches to doing 
this. One possibility is to evaluate \d a u n \^) by summa- 
tion over states as 

\daUnk) = |«mk) = • (51) 

This is very practical in the context of tight-binding cal- 
culations, where the sum over conduction bands runs 
only over a small number of terms, and we adopted this 
for the test-case calculations reported below. However, 
in first-principles calculations the sums over conduction 
states can be quite tedious, and one has to be careful 
to use the correct form for the velocity operator in the 
matrix elements (see discussion following Eq. (5)). Al- 
ternatively, the needed derivatives of |u n k) can be ob- 
tained from finite difference methods by making use of 
the discretized covariant derivative 33,34 as discussed in 
Appendix A. It may also be possible to use standard 
linear-response methods 35 to compute \d a u n -k), as this 
is an operation which is already implemented as part of 
computing the electric-field response in several modern 
code packages. 
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FIG. 2: 2x2 four-site square lattice used in the numerical 
tests. The absolute value of the hopping parameter t is set to 
1. $i. .4 are the threading fluxes through the four plaquettes. 



not all flux patterns break TR symmetry. For instance, 
for the flux patterns (<&i, $2, $3, $4) = {+</>, +<j>, —<t>, -<t>) 
and (+<j), —(f), +<j), —(f>), TR symmetry is restored by some 
spatial symmetry (an additional translational symmetry 
and a mirror symmetry, respectively); the orbital mag- 
netization then vanishes for any value of the parameter 
<j>. On the other hand, the flux pattern (2(p,—cj),0,—^)) 
violates inversion and mirror symmetry, and therefore re- 
alizes TR symmetry breaking. 

The on-site energies (Ea, Eb, Ec, Ed) have been set 
to the values (—3, 0, —3, 0). This choice results in an insu- 
lator with two groups of two entangled bands as shown in 
Fig. 3. Switching on the fluxes splits the bands along the 
X-L line, which are otherwise two-fold degenerate. The 
k-derivative of Bloch orbitals was computed by the sum- 
over-states formula Eq. (51). We treated the two lowest 
bands as filled and we verified that the multi-band Chern 
number is zero. 



A. Normal insulating case 

We present in this section numerical tests using a 
nearest- neighbor tight- binding Hamiltonian on a 2x2 
square lattice in which the primitive cell comprises four 
plaquettes, as shown in Fig. 2. This results in a four- 
band model. The modulus t of the (complex) nearest- 
neighbor hopping amplitude is set to 1, thus fixing the 
energy scale. TR breaking is achieved by endowing some 
of the hopping amplitudes with a complex phase fac- 
tor e 1 ^. This amounts to threading a pattern of mag- 
netic fluxes through the interiors of the four plaquettes, 
as shown in Fig. 2, in such a way that the threading 
flux $j is just the sum of the phase factors associated 
with the four bonds delineating plaquette i, counted 
with positive signs for counterclockwise-pointing bonds 
and minus signs for clockwise ones. The constraint 
of vanishing macroscopic magnetic field corresponds to 
$1 + $2 + $3 + $4 = 2-7T x integer. We found that 



<|> = jc/10 




FIG. 3: Band structure of the square lattice for <j> = 7r/10. 
The flux pattern is ($1, $2, $3, $4) = (2<f>, —<f>, 0, —<j>), and the 
on-site energies are (Ea,Eb,Ec,Ed) = (—3,0,-3,0) (see 
also Fig. 2). The two lower bands are treated as occupied. 
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FIG. 4: Orbital magnetization of the square-lattice model 
as a function of the parameter </>. The two lower bands are 
treated as occupied. Open circles: extrapolation from finite- 
size samples. Solid line: discretized k-space formula, Eq. (31). 



We built square finite samples, cut from the bulk, made 
of L x L four-site unit cells and having 2L + 1 sites on 
each edge. Their orbital magnetization (dipole per unit 
area) M(L) is straightforwardly computed as in Eq. (4). 
We expect the L — > oo asymptotic behavior 



M(L) = M + a/L + b/L 2 



(52) 



where M is the bulk magnetization according to Eq. (48). 
The terms a/L and b/L 2 account for edge and corner 
corrections, respectively. 

We performed calculations up to L = 14 (841 lattice 
sites). The resulting orbital magnetization as a function 
of the parameter (f> is shown in Fig. 4. We independently 
computed the bulk orbital magnetization M from a dis- 
cretization of the reciprocal-space formula Eq. (48). We 
get well converged results (to within 0.1%) for a 50x50 
k-point mesh in the full BZ. 

So far, we have studied a model multi-band insula- 
tor, having zero Chern number. For this specific case we 
provided above a solid analytic proof of our reciprocal- 
space formula, which holds in the thermodynamic limit. 
Indeed, the numerical results confirm the correctness of 
the k-space formula, while also providing some informa- 
tion about actual finite-size effects and numerical conver- 
gence. 



order to smooth Fermi-surface singularities, and to ob- 
tain well converged results, we adopt the simple Fcrmi- 
Dirac smearing technique, widely used in first-principle 
electronic-structure calculations. This amounts to re- 
place, the (integer) Fermi occupation factor 0(/z — e n k) 
with a suitable smooth function / M (e„k). We therefore 
replace in Eq. (48): 



2 -^S/M^nk)- 



(53) 



Reasoning in terms of a fictitious temperature, one may 
choose a Fermi-Dirac distribution 



1 



1 + exp[(e - fj,)/a] 



(54) 



In all subsequent calculations, we set a = 0.05 a.u., which 
provides good convergence. 

We compute the orbital magnetization as a function of 
the chemical potential /z with <f) fixed at ir/3. Using the 
same procedure as in the previous section, we compute 
the orbital magnetization by the means of the heuris- 
tic k-space formula, Eq. (48), and we compare it to the 
extrapolated value from finite samples, from L=8 (289 
sites) to L=16 (1089 sites). We verified that a k-point 
mesh of 100 x 100 gives well converged results for the bulk 
formula, Eq. (48). 

The orbital magnetization as a function of the chemical 
potential for cj> = tt/3 is shown in Fig. 5. The resulting 
values agree to a good level, and provide solid numerical 
evidence in favor of Eq. (48), whose analytical proof is 
still lacking. The orbital magnetization initially increases 
as the filling of the lowest band increases, and rises to a 
maximum at a fi value of about —4.1. Then, as the filling 
increases, the first (lowest) band crosses the second band 
and the orbital magnetization decreases, meaning that 
the two bands carry opposite-circulating currents giving 
rise to opposite contributions to the orbital magnetiza- 
tion. The orbital magnetization remains constant when 
fi is scanned through the insulating gap. Upon further 
increase of the chemical potential, the orbital magneti- 
zation shows a symmetrical behavior as a function of ji, 
the two upper bands having equal but opposite dispersion 
with respect to the two lowest bands (see Fig. 3). 



C. Chern insulating case 



B. Metallic case 

In the previous section we addressed the case of a TR- 
broken multi-band insulator, by treating the two low- 
est bands as occupied. Here we are going to extend our 
analysis to the metallic case. We are using the same 
model Hamiltonian as in the previous section, but we 
allow the Fermi level to span the energy range roughly 
from —5.45 to 2.45 energy units, namely from the bot- 
tom of the lowest band to the top of the highest one. In 



In order to check the validity of our heuristic Eq. (48) 
for a Chern insulator, we switch to the Haldane model 
Hamiltonian 11 that we used in a previous paper 7 to ad- 
dress the C = insulating case. In fact, depending on 
the parameter choice, the Chern number C within the 
model can be either zero or nonzero (actually, ±1). 

The Haldane model is comprised of a honeycomb lat- 
tice with two tight-binding sites per cell with site energies 
±A, real first-neighbor hoppings t\, and complex second- 
neighbor hoppings tie^ 1 ^ , as shown in Fig. 6. The result- 
ing Hamiltonian breaks TR symmetry and was proposed 
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FIG. 5: Orbital magnetization of the square-lattice model as a 
function of the chemical potential /i for <f> = 7r/3. The shaded 
areas correspond to the two groups of bands. Open circles: 
extrapolation from finite-size samples. Solid line: discretized 
k-space formula, Eq. (48). 
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FIG. 7: Orbital magnetization of the Haldane model as a 
function of the parameter <j>. The lowest band is treated as 
occupied. Open circles: extrapolation from finite size sam- 
ples. Solid line: Eq. (48). The system has non-zero Chern 
number in the region in between the two vertical lines. 
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FIG. 6: Four unit cells of the Haldane model. Filled (open) 
circles denote sites with Eq = — A (+A). Solid lines connect- 
ing nearest neighbors indicate a real hopping amplitude ti; 
dashed arrows pointing to a second-neighbor site indicates a 
complex hopping amplitude tie!''* '. Arrows indicate sign of the 
phase (j> for second-neighbor hopping. 



(for C = ±1) as a realization of the quantum Hall effect 
in the absence of a macroscopic magnetic field. Within 
this two-band model, one deals with insulators by taking 
the lowest band as occupied. 

In our previous paper 7 we restricted ourselves to C = 
to demonstrate the validity of Eq. (48), which was also 
analytically proved. In the present work we address the 
C =/= insulating case, where instead we have no proof of 
Eq. (48) yet. We are thus performing computer experi- 
ments in order to explore uncharted territory. 

Following the notation of Rcf. 11, we choose the pa- 
rameters A = 1, ti = 1 and \t 2 \ = 1/3. As a function of 
the flux parameter <j>, this system undergoes a transition 
from zero Chern number to \C\ — 1 when | sin<^| > l/\/3. 

First we checked the validity of Eq. (48) in the Chern 
insulating case by treating the lowest band as occupied. 
We computed the orbital magnetization as a function of 
4> by Eq. (48) at a fixed /i value, and we compared it to 
the magnetization of finite samples cut from the bulk. 



For the periodic system, we fix ti in the middle of the 
gap; for consistency, the finite-size calculations are per- 
formed at the same /i value, using the Fermi-Dirac dis- 
tribution of Eq. (54). The finite systems have therefore 
fractional orbital occupancy and a noninteger number of 
electrons. The biggest sample size was made up of 20 x 20 
unit cells (800 sites). The comparison between the finite- 
size extrapolations and the discretized k-space formula is 
displayed in Fig. 7. This hcuristically demonstrates the 
validity of our main results, Eqs. (46) and (48), in the 
Chern-insulating case. 

Next, we checked the validity of Eq. (48) for the most 
general case, following the transition from the metallic 
phase to the Chern insulating phase as a function of the 
chemical potential /x. To this aim we keep the model 
Hamiltonian fixed, choosing <j> = Q.7n; for /i in the gap 
this yields a Chern insulator. The behavior of the mag- 
netization while fi varies from the lowest-band region, to 
the gap region, and then to the highest-band region is 
displayed in Fig. 8, as obtained from both the finite-size 
extrapolations and the discretized k-space formula. This 
shows once more the validity of our heuristic formula. 
Also notice that in the gap region the magnetization is 
perfectly linear in /i, the slope being determined by the 
lowest-band Chern number according to Eq. (49). 



IV. CONCLUSIONS 

We present here a formalism for the calculation of the 
orbital magnetization in extended systems with broken 
TR symmetry, in the case of vanishing (or commensu- 
rate) macroscopic B field. This extends our previous 
work of Ref. 7 to the multi-band case, essential for real- 
istic calculations. 

First, we consider the case of zero Chern invariant, 
where we provide an an analytic proof, based upon the 
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Wannier representation. Our main result, Eq. (31), takes 
the form of a BZ integral of a gauge-invariant quan- 
tity, which can easily be computed using reciprocal- 
space discretization. We provide numerical tests for a 
two-dimensional model, where our discretized formula is 
checked against calculations performed for finite sam- 
ples cut from the bulk, with "open" boundary condi- 
tions. Our numerical tests appear to confirm that indeed 
Eq. (31) is the correct expression for the orbital magne- 
tization in a periodic system. 

Second, we propose a heuristic extension of Eq. (31) to 
the case of non-zero Chern invariant, based on the obser- 
vation that the integrand in Eq. (31) is gauge invariant, 
contrary to the analogous electrical case, where only the 
BZ integral is gauge- invariant, not the integrand. 8,9 On 
the basis of general considerations (namely, invariance by 
translation of the energy zero), the minimal modification 
extending Eq. (31) to the nonzero-Chern-number case 
yields Eq. (46). Remarkably, Eq. (46) is essentially iden- 
tical to a recent expression derived by Xiao et al. 10 in the 
context of a scmiclassical approximation. We check the 
full quantum-mechanical validity of Eq. (46) on a two- 
dimensional model by means of numerical tests, compar- 
ing to finite size calculations as above. The agreement 
is excellent, thus providing strong support for our for- 
mula, well beyond the semiclassical regime, even though 
we cannot yet provide an analytic proof of it. 

Third, since our heuristic Eq. (46) is well-defined for 
a KS metal, we also check the validity of Eq. (46) using 
the same two-dimensional model as for the metallic case, 
this time allowing the chemical potential /i to be varied 
through the bands. Once more the agreement is excel- 
lent, providing a numerical demonstration of the validity 
of Eq. (46). 

The electrical analogue of the present theory is the 
modern theory of polarization, 8 ' 9 developed in the 1990s, 
and valid for insulators only. When comparing that the- 
ory with the present one, in the insulating case, there is 
an important difference which is worth stressing. In the 
electrical case, the whole electronic contribution to the 
macroscopic polarization can be expressed in terms of the 
electric dipoles of the bulk WFs. This has a precise coun- 
terpart here, where the local-circulation contribution can 
in fact be expressed in terms of the magnetic dipoles of 
the bulk WFs. However, we have shown that in the mag- 
netic case there is an additional "itinerant-circulation" 
contribution which has no electrical analogue. When an- 
alyzing finite samples, the latter contribution appears to 
be due to chiral currents circulating at the sample bound- 
aries. Nonetheless, one of our major findings is that even 
this contribution can be expressed as a bulk, boundary- 
insensitive term. 

Both our original expression, Eq. (31), and its heuris- 
tic extension, Eq. (46), for the orbital magnetization of 
a crystalline solid can easily be implemented in existing 
first-principle electronic structure codes, making avail- 
able the computation of the orbital magnetization in 
crystals, at surfaces and in reduced dimensionality solids 




-0.6 1 ' ' ' ' ' ' ' ' 1 
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FIG. 8: Orbital magnetization of the Haldane model as a 
function of the chemical potential /i for cj> = 0.7-7T. The shaded 
areas correspond the position of the two bands. Open circles: 
extrapolation from finite-size samples. Solid line: discretized 
k-space formula, Eq. (48). 



such as nanowires. 
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Appendix A: Finite difference evaluation of the 
Chern invariant and magnetization 

Using the definition of the covariant derivative 33 ' 34 

\d a u n k) = Q k |<9 Q ?i„k), (55) 
Eqs. (33-35) can be rewritten as 

fa.af} = ^(5U"nk|<9/3U n k), (56) 

n 

gk,af3 = y^jdgUnklHuldpUnk), (57) 

n 

hk,aP = y^Tin'k (<9aU„'k|<9/3U nk ). (58) 
nn' 

We assume that the occupied wavefunctions |u„k) have 
been computed on a regular mesh of k-points, and we let 
bi, \>2, and t)3 be the primitive reciprocal vectors that 
define the k-mesh. Then the covariant derivative in mesh 
direction i can be defined as 

l^iWnk) = b iQ |d a u„ k ) (59) 
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(sum over a implied). Inserting this into Eqs. (56-58) and 
taking the antisymmetric imaginary part as in Eq. (41), 
we obtain 

f k = - eiji hi ^2 lm{djUnk\diu n k) , (60) 

n 

gk = - tijl b 4 ^ lm (djUnk\Hk\dlU n k) , (61) 

n 

hk = ^ eiji h t ^2 £ ; nr l 'kIni(9,M„'k|9/u nk ), (62) 

nn' 

where a sum over ijk is implied and v is the volume of 
the unit cell of the k-space mesh. On this mesh, the BZ 
integral in Eq. (42) becomes a summation 

C = 2^ X! e *J* bi Im (9jU„ k |9;u„ k ) (63) 

k n 

and similarly for the magnetization in Eq. (43). 

The appropriate finite-difference discretization of the 
covariant derivative in mesh direction i is 33 ' 34 

\diU n k) = 2 ( l"«,k+b t ) - |w„,k-bi)) (64) 

where \u n ,k+q) is the "dual" state, constructed as a lin- 
ear combination of the occupied |M nj k+q) at neighboring 
mesh point q, having the property that (w n 'k|w nj k+q) = 
8 n i n . This ensures that {u n 'k\diU n k) — consistent with 
Eq. (55), and is solved by the construction 33,34 

\Un,k+d) = (^k.k+q)"'" K',k+q) (65) 
n' 

where 

(5'k,k+q)rm' = ("tlk | U n > ,k+q) • (66) 

Eqs. (60-66) provide the formulas needed to calculate 
the three gauge-invariant quantities fk, gk, and hk 
on each point of the k-mesh. By summing these as 
in Eq. (63) it is straightforward to obtain C, Mlc, 
and Mic, respectively. Since we have derived this 
finite-difference representation using gauge-invariant 
quantities at each step, it is not surprising that we 



obtain the gauge- invariant contributions Mlc and Mic, 
as opposed to the gauge-dependent Mlc and Mic, from 
this approach. 



Appendix B: The nonAbelian Berry curvature 

It has been noticed in Sec. II D that the vector quantity 
fk is the Berry curvature. From Eqs. (38) and (41), this 
can be regarded as the trace of the N b x N h matrix F k 
having vector elements 

Fk,nn> = i (<9 k W„k| X l^k^n'k) 

- »y^(dkM«kK»k) x (u m k\dkU n 'k) ■ (67) 

m 

This quantity is known within the theory of the geo- 
metric phase as the nonAbelian Berry curvature, 36 and 
characterizes the evolution of an iV^-dimensional mani- 
fold (here, the states |w„k)) hi a parameter space (here, 
k-space). The nonAbelian curvature is gauge-covariant, 
meaning that if the states are unitarily transformed as 

\u n k) -» ^2u nn >(k)\u n , k ), (68) 

n' 

then the matrix Fk transforms as 

F k ,„n' -> Ul n (k)Fk, mm >U m >n>(k). (69) 

mm' 

This implies that the invariants of the matrix Fk, such 
as its trace fk 2 are gauge- invariant. In fact, as discussed 
in Sec. II D, fk behaves like a standard (i.e., Abelian) 
curvature. 

We also notice that the energy matrix Ek, Eq. (1), is 
also gauge-covariant in the sense of Eq. (69). It is then 
easy to verify that the trace (over the band indices) of 
the matrix product Ek Fk is a gauge-invariant quantity. 
In fact, this trace is identical to hk as defined in Sec. II D, 
whose gauge-invariance we proved in a different way. The 
special Nb = 1 case was previously dealt with in Ref. 7, 
where the analogue of hk takes the form of the product 
of energy times curvature, both gauge-invariant quanti- 
ties. The present finding shows that, in the multi-band 
case, this must be generalized as the trace of the (matrix) 
product Ek times Fk, both gauge-covariant quantities. 
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